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Abstract: The objectives of the study are to integrate the conditional Latin Hypercube 
Sampling (cLHS), sequential Gaussian simulation (SGS) and spatial analysis in remotely 
sensed images, to monitor the effects of large chronological disturbances on spatial 
characteristics of landscape changes including spatial heterogeneity and variability. The 
multiple NDVI images demonstrate that spatial patterns of disturbed landscapes were 
successfully delineated by spatial analysis such as variogram, MoranT and landscape 
metrics in the study area. The hybrid method delineates the spatial patterns and spatial 
variability of landscapes caused by these large disturbances. The cLHS approach is applied 
to select samples from Normalized Difference Vegetation Index (NDVI) images from 
SPOT HRV images in the Chenyulan watershed of Taiwan, and then SGS with sufficient 
samples is used to generate maps of NDVI images. In final, the NDVI simulated maps are 
verified using indexes such as the correlation coefficient and mean absolute error (MAE). 
Therefore, the statistics and spatial structures of multiple NDVI images present a very 
robust behavior, which advocates the use of the index for the quantification of the 
landscape spatial patterns and land cover change. In addition, the results transferred by 
Open Geospatial techniques can be accessed from web-based and end-user applications of 
the watershed management. 
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1. Introduction 

The interest in land cover analysis at regional to global scales has grown dramatically in the last 
decade [1]. To ensure a sustainable management of natural resources, it is necessary to understand and 
quantify the landscape change processes. Patterns of landscape change are the results of complex 
interactions between physical, biological and social forces [2]. To understand and predict change 
processes, one needs to monitor and characterize spatial patterns of land use/land cover 
change [3]. Therefore, the methods for proving the sufficient information to understand the land cover 
change, such as sampling and simulation are very important. Remote sensing has emerged as the most 
useful data source for quantitatively measuring land-cover changes at the landscape scale. The 
dynamics of change processes can be investigated through temporal series of remote sensing data and 
by analyzing change trajectories [3,4]. Remotely sensed data can describe surface processes, including 
landscape dynamics and land cover change, as such data provide frequent spatial estimates of key 
environmental variables [5]. Land cover change is regarded as the single most important variable of 
global change affecting ecological systems with an impact on the environment [6]. Satellite image data 
have become an important source of information for monitoring the vegetation conditions, land use 
and land cover change [7]. The most widely used vegetation index in this context is NDVI, the 
normalized difference vegetation index, which is a function of red and near-infrared spectral bands [8]. 
The NDVI value indicates a level of photosynthetic activity [9]. The change areas can be identified 
through the subtraction of the NDVI image of one date from the NDVI image of another date [10]. 
Significant differences of NDVI images in landscape changes, such as landslides are induced by a 
disturbance. In addition, the landscapes are a composite of complex nature and biophysical forces that 
are manifested on the landscape through the composition and spatial organization of land cover and 
variations in NDVI. Multi-temporal NDVI images are practical for monitoring vegetation dynamics on 
a regional scale. 

The land cover change areas usually cover only a small portion of the entire image and may not be 
detected or represented by simple random sampling or systematic sampling unless the sampling 
intensity is very high .The principle of sampling for land cover analysis is well established [1,10,11] . 
Random sampling is statistically appealing because of the applicability of the classical statistical 
procedures. But it is not generally an efficient approach. For this reason, other sampling designs such 
as stratified random, systematic grid, and others may be preferred for land cover analysis. However, 
the sampling problem is complicated by the nature of fine resolution satellite data. For reasons of costs 
and efficiency of using the acquired data, it is much preferable to select scenes that will make the 
greatest contribution to the characterization of land cover over the entire domain [1]. Thus, Latin 
hypercube sampling (LHS) was used to obtain reference data for this study. LHS was initially 
developed for the purpose of Monte-Carlo simulation, efficiently selecting input variables for the 
models [12,13]. Conditional Latin hypercube sampling (cLHS), which is based on the empirical 
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distribution of original data, provides full coverage of each variable by maximally stratifying the 
marginal distribution and ensuring a good spread of sampling points [14]. The cLHS is an effective 
way to replicate the distribution of variables compared with random sampling and stratified spatial 
sampling, and ensures the correlation of sampled variables to replicate the original data [14,15]. After 
cLHS sampling, geostatistical conditional simulations with the sampling data can be applied to 
simulate the spatial variability and spatial distribution of interest [16]. For sampled data, a 
geostatistical conditional simulation technique, such as sequential Gaussian simulation (SGS), can be 
applied to generate multiple realizations, including an error component, which is absent from classical 
interpolation approaches [17-20]. Sequential simulation procedures generate simulated values using a 
Monte Carlo approach at each node of the simulation grid visited sequentially according to a random 
path. SGS realization varies greatly in space on the condition of original data and all previously 
simulated values [19,21]. Moreover, development of efficient procedures for designing information- 
effective samplings and simulations is an essential task for more accurately understanding the 
spatial distribution. 

Spatial analysis (i.e. MoranT, and variogram) delineates spatial variations and patterns of land 
cover from remote sensing images [6,22,23]. Variograms are crucial to geostatistics. Both are 
functions related to the variance of spatial separation and provide a concise description of the scale and 
pattern of spatial variability [24]. Samples of remotely sensed data (e.g., satellite or air-borne sensor 
imagery) can be employed to construct variograms for remotely sensed research [24]. Moreover, 
MoranT and variograms have been used widely to understand the nature and causes of spatial variation 
within an image [22,25,26]. The MoranT and variogram characterize and quantify spatial 
autocorrelation, heterogeneous spatial components (spatial variability and spatial structure) of a 
landscape and the spatial heterogeneity of vegetation cover at the landscape level [23,25,26]. 
Furthermore, landscape metrics have been used increasingly to assess land-use change in the last 
decade [6,27-30]. Variables that characterize landscape patterns, such as the number, area, and spatial 
patterns of different patch types, change when the land use/ land cover change is altered. Tandscape 
composition, configuration, and connectivity are primary descriptors of landscape patterns [27], and 
can be quantified using spatial landscape indices or metrics. Tandscape composition refers to the 
characteristics associated with the variety and abundance of patch types within a given landscape. 
Spatial configuration of a landscape denotes the spatial characteristics and arrangement, position, or 
orientation of patches within a class or landscape [31]. In addition, landscape metrics have proved 
effective in land-use planning and design because they can characterize the landscape pattern in the 
design alternatives [6,32,33]. Accordingly, monitoring, delineating and sampling landscape changes, 
spatial structure and spatial variation induced by large physical disturbances are essential to landscape 
management and restoration, and disaster management in Taiwan. Earthquakes have long been 
recognized as a major cause of landslides [34,35]. The earthquake seriously damaged the vegetation 
and landslides in the region [35-37]. Additionally, typhoons are the main external triggering factor to 
debris flow and are extremely important natural disturbances that characterize the structure, function 
and dynamics of many ecosystems [3 8-40]. Typhoons that bring tremendous amounts of rainfall hit 
Taiwan every year from July to October [41]. During 1996-2004, large disturbances in the following 
sequence impacted central Taiwan: such as typhoon Herb (August 1996); the Chi-Chi earthquake 
(September 1999); typhoon Xangsane (November 2000); typhoon Toraji (July 2001); typhoon Dujuan 
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(September, 2003), and typhoon Mindulle (June 2004) [41]. After the Chi-Chi earthquake, the 
expansion rate of landslide areas increased 20-fold in central Taiwan [42]. The influences of large 
physical disturbances on ecosystem structure and function have garnered considerable attention [43-46]. 
Numerous extension cracks, which accelerate landslides during downpours, were generated on hill 
slopes during the Chi-Chi earthquake [37,40,47]. The landslides are the natural hazards and the 
resulting in a variety of human and eco-environmental impacts [20]. 

In the study, the NDVI data derived from SPOT images before and after the Chi-Chi earthquake in 
the Chenyulan basin of Taiwan, as well as images before and after four large typhoons (Xangsane, 
Toraji, Dujuan and Mindulle) were analyzed to identify the landscapes change caused by these major 
disturbances. For the representation of the changes of Multi-temporal NDVI data, we used the hybrid 
methods including spatial sampling, simulation, spatial analysis and landscape metrics in the 
space-time data analysis. The spatial analyses identify the landscapes change and delineate spatial 
variations of NDVI images before and after large physical disturbances in central Taiwan. Moreover, 
conditional LHS (cLHS) schemes with NDVI images were used to select spatial samples from actual 
NDVI images to detect landscape changes induced by a series of large disturbances. Furthermore, the 
best cLHS samples selected with the NDVI values were used to simulate NDVI distributions using 
SGS. Finally, the analyses are used to evaluate the approaches for defining a sufficient sampling size 
to capture the spatial variability of the land cover. The study can be used to understand the NDVI 
spatial structures within an image domain for the watershed management. 

2. Methods and Materials 

2.1. Study Area and Remote Sensing Data 

The Chenyulan watershed, located in central Taiwan, is a classical intermountain watershed which 
is traversed by the Chenyulan stream in the south to north direction. The average altitude and area of 
this watershed are 1,540 m and 449 km , respectively (Figure 1). Differences in uplifting along the 
fault generated abundant fractures over the watershed and resulted in an average slope of 62.5% and 
relief of 585 m/km 2 . Moreover, the main course of the Chenyulan stream had a gradient of 6.1%, and 
more than 60% of its tributaries had gradients exceeding 20%. In this area, slates and meta-sandstones 
are the dominant lithologies in the metamorphic terrains [25,42]. Based on the relative amounts of slate 
and meta-sandstone, the metamorphic strata in the eastern part of the study area are divided into four 
parts: Shihpachuangchi, Tachien Meta- Sandstone, Paileng Meta- Sandstone, and Shuichangliu [25,42]. 
Figure 1 also shows the six land cover categories in the area on November 19, 2004 which include 
forest, grassland, cultivated land, water, bare land, and built-up. The details could be referred to [48]. 
The study area with dimensions of 5 x 5 km 2 (250 x 250 pixels) was selected from the upstream of the 
large debris flood announced in the watershed, as shown in Figure 1 . 

In this paper, the NDVI changes in land cover reflect the natural consequences of typhoons and 
earthquakes. At 01:4712.6" on September 21, 1999, an earthquake of magnitude 7.3 on the Richter 
scale and with a focal depth of 8.0 km struck central Taiwan. The cause of this event, known as the 
Chi-Chi earthquake, was movement of the CheLungPu fault. The epicenter was located at 23.87°N and 
120.75°E, near the Chenyulan watershed in south Nantou County. The magnitude of the earthquake 
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was estimated to be ML = 7.3 (ML: Local Magnitude or Richter Magnitude), and the rupture zone, 
defined by the aftershocks, measured about 80 km north-south by 25-30 km downdip [47,49]. 
Iso-contour maps of the earthquake's magnitude were reproduced from the Central Weather Bureau 
(Figure 1) [50]. After the earthquake, the center of typhoon Xiangsane moved from south to north 
through eastern Taiwan from October 31, 2000 to November 1, 2000 [51], with a maximum wind 
speed of 138.9 km/hr and a radius of 250 km (Figure 1). The maximum daily rainfall was 550 mm/day. 
The Toraji typhoon swept across central Taiwan from east to west on July 30, 2001 [52] (Figure 1). 
The typhoon brought extremely heavy rainfall, from 230 to 650 mm/ day, and triggered numerous 
landslides in Taiwan. Typhoon Toraji became a tropical storm after crossing Taiwan; however it 
brought 339 to 757 mm of total accumulated rainfall in the watershed [52] (Figure 1). From August 31, 
2003 to September 2, 2003, typhoon Dujuan was unusual in the double eye with a maximum wind 
speed of 165.0 km/hr, a radius of 200 km and maximum rainfall 200 mm/hr. From June 29, 2004 to 
July 2, 2004, Mindulle with maximum wind speed of 200.0 km/hr, a radius of 200 km. Severe rainfall 
(with a maximum amount of 787 mm) occurred over central-southwestern Taiwan (Figure 1) [53]. 

Figure 1. Typhoon paths and Chi-Chi earthquake center are around Taiwan, and land use 
patterns of the study area. 



Grassland 
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2.2. NDVI 



The NDVI is an empirically derived index used to estimate plant biomass through the integration of 
the red-visible and near-infrared spectral regions to represent plant pigmentation and chlorophyll 
content, respectively, in the characterization of land cover conditions [28]. Seven cloud- free SPOT 
images (1996/11/08, 1999/03/06, 1999/10/31, 2000/11/27, 2001/11/20, 2003/12/17 and 2004/11/19) of 
the Chenyulan watershed were purchased from the Space and Remote- sensing Research Center, 
Taiwan. The NDVI images of the study area were generated from SPOT HRV images with a 
resolution of 20 m according to the following equation: 

NIR + R 

where NIR and R are near-infrared and visible-red spectral data, respectively. The NDVI values range 
from -1 to +1; a high NDVI value represents a large amount of high photosynthesizing vegetation [53]. 



2.3. Landscape Matrices 



A wide variety of indices developed to characterize landscapes has often been applied to study 
spatial landscape patterns in landscape ecology; they can be categorized into: area/density/edge, shape, 
core area, isolation/proximity, contrast, contagion, connectivity, and diversity metrics [31]. Since 
landscape metrics can be used to analyze the spatial and temporal changes in landscape patterns, they 
provide a quantitative approach for studying the land cover change through the measurement of spatial 
and temporal variations in these metrics. To assess changes in land-use patterns, the following 
landscape metrics were calculated using the Patch Analyst in the GIS software Arc View 3.2a [55] is 
designed to compute a wide variety of landscape metrics for categorical map patterns. In this study, the 
following metrics were used (Table 1): the number of patches (NP), mean patch size (MPS), Patch 
Size Standard Deviation (PSSD), Patch Size Coefficient of Variance (PScov), Edge Density (ED), 
mean shape index (MSI) and mean nearest neighbor (MNN) [31]. 



2.4. Variogram 



In geostatistical methods, variograms can be used to quantify the observed relationship between the 
values of samples and the proximity of samples [19]. Following the work of Garrigues et al. [23], 
Garrigues et al. [26] and Lin et al. [35], NDVI data are considered values of punctual regionalized 
variable. An experimental variogram for interval lag distance class h, y{h) , is represented by: 



7 (h) = -—y j [Z(x i+ h)-Z(x i )f (2) 
2n(h) , =1 v > 

where h is the lag distance that separates pairs of points; Z(x) is NDVI value at location x, and Z(x + h) 
is NDVI value at location x + h; n(h) is the number of pairs separated by lag distance h. 
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Table 1. Landscape metrics. 
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Edge Density (ED) 
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Mean shape index (MSI) 
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Shape metrics 












Mean nearest neighbor 
(MNN) 


MNN = 7=1 






Diversity 
metrics 



where rii is the number of patches in land-use class i; ay is the j patch area (ha.) inland-use class i; 
m is the total number of patch classes; e ik is the total length (m) of the edge between patch classes / 
and k\ py is the f h patch perimeter (m) in land-use class i; hy is the distance (m) from the / h patch to 
the nearest neighboring patch of the same class /, based on the edge-to-edge distance. 



2.5. Conditional Latin hyper cube 

The cLHS procedure represents an optimization problem: given N sites with ancillary data (Z), 
select n sample sites (n « N) such that the sampled sites form a Latin hypercube. For k continuous 
variables, each component of Z (size, N x k) is divided into n (sample size) equally probable strata 
based on their distributions, and z (size n x k) is a sub-sample of Z. The procedures of the cLHS 
algorithm [14] are as follows: 

1 . Divide the quantile distribution of Z into n strata, and calculate the quantile distribution for each 
variable, q l j9 ... 9 q" +l . Calculate the correlation matrix of Z (C). 

2. Pick n random samples from N; calculate the correlation matrix of z (T). 
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3. Calculate the objective function. The overall objective function combines to three components of 
the objective function (0\ 9 02, and O3). For general applications, all weightings are set to equal 
for all components of the objective function. 

A. For continuous variables, 

where r/(q l j < z. < q l ; +l ) is the number of zj that falls between quantiles q l . and q 1 * 1 

B. For categorical data, the objective function is to match the probability distribution for each 
class of 



-k. 



J 



(4) 



where r\ '(zj) is the number of z that belongs to class j in sampled data, and kj is the proportion of 
class j in Z. 

C. To ensure that the correlation of the sampled variables will replicate the original data, another 
objective function is added: 

°3=i;i>,-'„i (5) 

i=\ j=\ 

where c is the element of C, the correlation matrix of Z, and t is the equivalent element of T, the 
correlation matrix of z. 

4. Perform an annealing schedule:M = exp [-A0/T], where AO is the change in the objective 
function, and T is a cooling temperature (between 0 and 1), which is decreased by a factor d 
during each iteration. 

5. Generate a uniform random number between 0 and 1. If rand < M, accept the new values; 
otherwise, discard changes. 

6. Try to perform changes: Generate a uniform random number rand. If rand < P 9 pick a sample 
randomly from z and swap it with a random site from unsampled sites r. Otherwise, remove the 
sample(s) from z that has the largest r/(q l j < z. < q 1 ^ 1 ) and replace it with a random site(s) from 

unsampled sites r. End when the value of P is between 0 and 1, indicating that the probability of 
the search is a random search or systematically replacing the samples that have the worst fit with 
the strata. 

7. Go to step 3. Repeat steps 3-7 until the objective function value falls beyond a given stop 
criterion or a specified number of iterations. 



2. 6. Sequential Gaussian Simulation 



The SGS assumes a Gaussian random field, such that the mean value and covariance completely 
characterize the conditional cumulative density function [56]. During the SGS process, Gaussian 
transformation of available measurements is simulated, such that each simulated value is conditional 
on original data and all previously simulated values [21,57]. A value simulated at a one location is 
randomly selected from the normal distribution function defined by the kriging mean and variance 
based on neighborhood values. Finally, simulated normal values are back-transformed into simulated 



Sensors 2009, 9 



6678 



values to yield the original variable. The simulated value at the new randomly visited point value 
depends on both original data and previously simulated values. This process is repeated until all points 
have been simulated. In sequential simulation algorithm, modeling of the N-point cumulative density 
function (ccdf) is a sequence of N univariate ccdfs at each node (grid cell) along a random path [58]. 
The sequential simulation algorithm has the following steps [58]: 

1. Establish a random path that is visited once and only once, all nodes {jc,-, / = 1, A, N} 
discretizing the domain of interest Doman. A random visiting sequence ensures that no spatial 
continuity artifact is introduced into the simulation by a specific path visiting N nodes. 

2. At the first visited N nodes x\ : 

A. Model, using either a parametric or nonparametric approach, the local ccdf of Z(x\) 
conditional on n original data {Z (x a ) 9 a=l 9 A, n} F z (x\\ z\\{n)) = prob {Z (x i) < z\\(n)} 

B. Generate, via the Monte Carlo drawing relation, a simulated value x \) from this ccdf 
Fz(x\ \ z\\(n)) 9 and add it to the conditioning data set, now of dimension n + 1, to be used 
for all subsequent local ccdf determinations. 

3. At the i t h node jc,- along the random path: 

A. Model the local ccdf of Z(x i) conditional on n original data and the / - 1 near previously 
simulated values { z (/) (x/),y = 1,A, / - 1}: 

F z (x x ; z. \{n + i - 1)) = prob{Z(x i ) < z. \(n + i - 1)} (6) 

B. Generate a simulated value z {l) {x l ) from this ccdf and add it to the conditioning data set, 
now of dimension n + /. 

4. Repeat step 3 until all N nodes along the random path are visited. 

2. 7. Moran 's I 

Spatial autocorrelation is a useful tool for describing the dependency of spatial patterns. First, 
spatial structures are described by so-called structure functions [25,59].Moran's I, which ranges 
between -1 and +1, is a well known spatial autocorrelation method [60]. The index, I, is calculated as 
follows: 

a/«)Z>,-^) 2 

where yh and yi denote the values of the observed variable at sites h and I, respectively; and Whi denotes 
the weight of the variable. The weights, wy, are written in an (n x n) weight matrix W, which is the 
sum of the weights w^ for a given distance class [61]. Moran' s I is high and positive when a value is 
similar to adjacent values, and low when a value is dissimilar to adjacent values. In this paper, the 
global Moran' s I value for the NDVI was calculated to compare the spatial relations of the NDVI 
among various events. As a result, the phenomenon of spatial autocorrelation of NDVI could be tested. 
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3. Results and Discussion 

3.1. Statistics and Spatial Analysis of NDVI Images 

The NDVI is one of the most popular methods for monitoring vegetation conditions. It has been 
reported that multitemporal NDVI is useful for classifying land cover and the dynamics of 
vegetation [19,62,63]. However, the typhoons and earthquakes is a major natural disturbance to land 
cover change in Taiwan. For example, the Chi-Chi earthquake led to landslides, dammed lakes and a 
high death toll. Like the typhoons, subsequent rainstorms cause divergent destruction of vegetation; 
this destruction may be influenced by the rainfall distribution and typhoon path. However, due to 
destruction of vegetation and increased soil exposure in the Chenyulan watershed, the landslide ratio 
increased with successive rainstorms and strong earthquakes and resulted in the decreased NDVI 
values [41]. 

Figure 2 shows the NDVI images of the area on seven events of the large disturbances. Table 2 and 
Figure 3 summarize the statistics and histograms for seven actual NDVI images of the area before and 
after disturbances. The lowest mean and minimum NDVI values in 1996-2004 occurred on October 31, 
1999, after the Chi-Chi earthquake in the studied area. Moreover, the largest range between minimum 
and maximum NDVI values also occurred on October 31, 1999, after the Chi-Chi earthquake in the 
study area. The other negative minimum NDVI values occurred on November 27, 2000, and December 
17, 2003, in the study area. The second and third greatest impacts on all landscapes are from typhoons 
Xangsane (November 2000) and Dujuan (September 2003) in the area, respectively (Figure 2 and 
Table 2). During the dates, the standard deviations of NDVI values were slightly larger than those on 
the other dates. These statistical results illustrate that the Chi-Chi earthquake had the largest impact on 
all landscapes represented by NDVI images for the study area. Unfortunately, typhoon Xangsane that 
came after the Chi-Chi earthquake was the second disturbance to impact landscape changes in the 
study area. Numerous extension cracks, which increase the number of landslides during downpours, 
were generated on hill slopes during the Chi-Chi earthquake [25]. Statistics of remotely sensed images 
can be used as a basic tool to characterize landscape changes [64-68]. Statistical results illustrate that 
the effects of disturbances on the watershed landscape in the study area were significantly different in 
space and time over the entire landscape. 



Table 2. Statistics of NDVI images. 



Date 


Mean 


Std. 


Min. 


Max. 


Skewness 


Kurtosis 


1996/11/08 


0.36 


0.04 


0.11 


0.48 


-0.98 


1.45 


1999/03/06 


0.32 


0.04 


0.13 


0.43 


-0.58 


0.08 


1999/10/31 


0.14 


0.07 


-0.22 


0.33 


-1.23 


1.35 


2000/11/27 


0.15 


0.07 


-0.14 


0.35 


-0.47 


-0.30 


2001/11/20 


0.37 


0.05 


0.03 


0.50 


-1.03 


1.34 


2003/12/17 


0.15 


0.06 


-0.12 


0.33 


-0.27 


0.00 


2004/11/19 


0.35 


0.06 


0.05 


0.54 


-0.44 


0.07 
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Figure 2. NDVI images of the area on (a) 1996/11/08, (b) 1999/03/06, (c) 1999/10/31, (d) 
2000/1 1/27, (e) 2001/1 1/20, (f) 2003/12/17, and (g) 2004/1 1/19. 
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Figure 3. Histogram of original data on (a) 1996/11/08, (b) 1999/03/06, (c) 1999/10/31, (d) 
2000/11/27, (e) 2001/11/20, (f) 2003/12/17, and (g) 2004/11/19. 
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The first three lowest mean NDVI images in 1996-2004 occurred on 1999/10/31 (Table 2), after the 
Chi-Chi earthquake, typhoons Xangsane (2000/11/27), and Dujuan (2003/12/17) in the study area. 
Figure 4 shows mean NDVI varying with time and the classification data in the events for the low 
mean NDVI values and the brown parts in the images represent sparse vegetation when the NDVI 
values are less than zero. Landscape metrics were used to assess the damage of the three events include 
the Chi-Chi earthquake, typhoon Xangsane and Dujuan, and demonstrate the landscape metrics of 
land-use patterns in the low NDVI class level. The results from mean NDVI calculation revealed a 
stable plant growth and vegetation recovery tendency of the study area. These results confirm the 
previous studies that the vegetation keeps on recovering in the landslide areas but earthquake and 
subsequent rainstorms may impact on the vegetation recovery rates [20,69,70]. Figure 5 shows 
landscape matrices for (a) the number of patches (NP), (b) mean patch size (MPS), (c) Patch Size 
Standard Deviation (PSSD), (d) Patch Size Coefficient of Variance (PScov), (e) Edge Density (ED), (f) 
Mean Shape Index (MSI), (g) Mean Nearest Neighbor (MNN) in the events. Form the analysis of the 
indexes, the NP and ED for low NDVI class of the earthquake is larger than that of typhoon events. It 
is discovered that there is a large number of patches and edge density in the earthquake. The value of 
MPS for low NDVI class is the largest on 2000/11/27, middle on 1999/10/31, and the smallest on 
2003/12/17. The shape index (MSI) and edge index (ED) present robust behaviors, which advocate the 
use of the indexes for the quantification of the overall irregularity of patch shapes and edge in low 
NDVI class. The MSI for low NDVI class on 2003/12/17 (after typhoon Dujuan) is the smallest and 
represents that the landslide patch shape is close to the square. Moreover, the MNN for low NDVI on 
1999/10/31 (after the Chi-Chi earthquake) is less than that of the typhoon events. Classification data of 
NDVI used within landscape metrics also have a number of advantages as them give us information on 
several indices of spatial distribution including patch size, edge, and shape of patches, all of which are 
theoretically capable of providing different insights on landscape fragmentation after the large 
chronological disturbances. These indices show the fragmented NDVI patterns due to the natural 
disturbances in the watershed. Accordingly, the earthquake induces the major and regional damage but 
the disturbances of NDVI caused by the typhoons are minor and local in the study area. 

Previous studies that quantified the impact of large disturbances did not evaluate the spatial 
structures of NDVI images in the study areas. To depict spatial autocorrelation, landscape 
heterogeneity, spatial variability and patterns, MoranT, experimental variograms and their variogram 
models were first analyzed and fitted to seven images of the studied area (Figure 6 and Table 3) [20]. 
Figure 6(a) shows correlograms of the Moran's I (distance as the upper distance of a lag) in the seven 
events. All show positive spatial autocorrelation in short distance and negative ones in large distance, 
which is called gradient spatial pattern. Typhoons and earthquakes will disturb the spatial 
autocorrelation of the NDVI. Spatial correlations for all events among watershed land patches at 
distances on less than 2700 m are considerably positive. In the figure, the spatial autocorrelation of 
NDVI images in the area are highest on 1999/03/06 and are the lowest on 1999/10/31. Between the 
two dates, the Chi-Chi earthquake happened. Hence, the earthquake will cause the most serious 
disturbance and reduce the spatial autocorrelation. In Figure 6(b), the variogram models of the seven 
NDVI images for the study area are exponential models. The Sill (Co + Q of NDVI images from high 
to low in the area are on 1999/10/31, 2000/11/27, 2004/11/19, 2003/12/17, 2001/11/20, 1999/03/06 
and 1996/11/08. The sill is the upper limit that a variogram approaches at a large distance, and is a 
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measure of the variability of the investigated variable: a higher sill corresponds to greater variability in 
the variable. The spatial variations of NDVI images increase considerably from 1996/11/08 to 
1999/10/31 (after the Chi-Chi earthquake) in the area. The Nugget effects of NDVI images on 
1999/10/31 are larger than those in other images. The nugget effect is exhibited by the apparent non- 
zero value of the variogram at the origin, which may be due to the small-scale variability of the 
investigated process and measured errors. Exponential models with large sills and large nugget effects 
NDVI images are indicative of significant spatial heterogeneous landscapes induced by the Chi-Chi 
earthquake in the area. Moreover, typhoons Xangsane and Dujuan also generated heterogeneous 
landscapes in the area [20]. Spatial correlations of NDVI over the watershed at distances less than 
2,800 m are considerably positive. Spatial autocorrelation (Moran's I) was used to delineate spatial 
variations in the landscape and class levels of the sub-watersheds before and after disturbances [25]. In 
the paper, the spatial correlation tendencies are similar but the changes before and after disturbances 
are significantly different. Moreover, variography results confirm that the impacts of disturbances on 
the watershed landscape pattern were cumulative, but were not always evident in space and time in the 
entire landscape [9,71]. Variography results illustrate that NDVI discontinuities between fields create a 
mosaic spatial structure resulting primarily from large disturbances, such as the Chi-Chi earthquake, in 
the study area. In addition, the landscape metrics examine the patterns of landscape fragmentation after 
the large disturbances on the watershed landscape. The landscape metrics also indicate that the 
disturbances and disturbance regime are characterized by a variety of attributes, including size, 
frequency, intensity, severity, and shape [27]. 

In this study, spatial analysis and modeling results indicate that large disturbances, such as the Chi- 
Chi earthquake, created extremely complex heterogeneous patterns across the landscape. Thus, a 
disturbance may affect some areas and disturbance severity often varies considerably within an 
affected area on the landscape level [25,45]. However, the earthquake and the typhoons impacted the 
fragmentation, shape, isolation, and interspersion of the patches. These results verify that disturbances 
create complex heterogeneous patterns across the landscape because the severity of the disturbances 
frequently varies considerably within the affected area. Disturbances are not the only destructive and 
restorative causes of modification of the geomorphic landscape or the structure and composition of the 
forest that occur between disturbances [72]. The spatial analysis results of NDVI images are sufficient 
to present landscape changes induced by disturbances, particularly via spatial structure, variability and 
spatial correlation. Previous studies [25,41,70] indicated that landslides in the Chenyulan watershed 
were impacted by the Chi-Chi earthquake; however, the effect of the earthquake decreased as the time 
after the Chi-Chi earthquake [70]. Landslides induced by earthquakes and typhoons have distinct 
spatial patterns [41]. Bare land and grassland have high potential to become landslides during the large 
disturbances [48]. The high-magnitude Chi-Chi earthquake created these landscape variations in space 
in the Chenyulan watershed [25]. Moreover, typhoons significantly influence NDVI variations via the 
flow of accumulated rainfall and wind gradients [19]. Thus, the modeling results also prove the 
previous studies indicating that earthquake and subsequent rainstorms may cause divergent destruction 
of vegetation, and then this destruction may be influenced by the precipitation distribution and 
typhoon path. 
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Figure 5. Landscape matrices in (a) the number of patches (NP), (b) mean patch size 
(MPS), (c) Patch Size Standard Deviation (PSSD), (d) Patch Size Coefficient of Variance 
(PScov), (e) Edge Density (ED), (f) Mean Shape Index (MSI), (g) Mean Nearest Neighbor 
(MNN) in the event 1. (1999/10/31), event 2 (2000/1 1/27) and event 3 (2003/12/17). 
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Figure 6. (a) Moran'I and (b) experimental variograms of NDVI images before and after 
disturbances in the area. 
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Table 3. Variogram models of NDVI images. 



Date 


Model 




Parameters 




RSS 


r 2 


CO (me/ks) 2 


+C(m2/k2) 2 








1 QQ6/1 1 /OR 


Exp. 


0.000138 


0.001326 


654 


1.61E-08 


0.953 


1QQQ/0V06 
i y y y I \J u I \J\j 


Exp. 


0.000712 


0.001814 


4620 


6.07E-08 


0.945 


1999/10/31 


Exp. 


0.000590 


0.004440 


564 


1.68E-07 


0.939 


2000/1 1/27 


Exp. 


0.000186 


0.004676 


2646 


2.47E-07 


0.952 


2001/11/20 


Exp. 


0.000121 


0.002429 


1281 


5.62E-08 


0.933 


2003/12/17 


Exp. 


0.000126 


0.003126 


2298 


1.57E-07 


0.949 


2004/11/19 


Exp. 


0.000116 


0.003832 


1680 


1.19E-07 


0.977 



Exp.: Exponential model; Sph.: Spherical model; CO: Nugget; C0+C: Sill; R: Range; RSS: 
Residual Sums of Squares 



Furthermore, the results from mean NDVI revealed vegetation recovery tendency and a stable plant 
growth cycle for the study area. The preview study shows the same results that the vegetation recovery 
rate reached more than half of (58.93%) original vegetation regeneration in the landslide areas over 
two years of monitoring and assessing [69]. The poor recovery locations were distributed mainly in 
mountain ridge, scoured slope base and acidic sulfate soil areas. The areas had lower recovery rates 
affected by the impacts of the slope and the property of the soil [69,73]. 

3.2. Simulation with Selected Samples for Multiple Images 

cLHS is a stratified random procedure that provides an efficient way of sampling variables from 
their multivariate distributions. It provides a full coverage of range each variable by ensuring a good 
spread of the sampling points [12]. In the previous study [20], experimental variograms of cLHS 
samples with their NDVI values were constructed using the same lag interval to compare the spatial 
structures of the actual NDVI images. Lin et al. [20] listed the statistics and variogram for different 
samples from multiple NDVI images with 62,500 grids using the cLHS approach. The statistics for 
these that select more than 3,000 samples indicate the statistics obtained by cLHS can capture statistics 
of all actual NDVI images. The statistical and variogram analyses of cLHS samples also illustrate that 
the cLHS approach can be applied to select samples and capture the spatial structures of multiple 
historically accurate NDVI images [20]. The distributions of selected samples confirm that samples 
selected using cLHS provide a good coverage of the study area and are well spread and partially 
clustered in the study area [12]. These samples can be used in further monitoring and simulation to 
determine the impacts of disturbances on study landscapes in the latter. 

In this study, SGS simulation was performed based on the different samples for seven NDVI images 
in the area. Table 4 shows the statistics of the SGS simulations using 100, 500, 1,000, 2,000, 3,000, 
5,000, 7,000, and 10,000 cLHS selected samples in the seven events. The comparison reveals that the 
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statistics of these SGS simulations are matched to the original data. Over 3,000 samples, the statistics 
of these simulations are close to that of the original data. Figure 7 shows the average of 1,000 NDVI 
realizations produced by SGS simulations with 100, 500, 1,000, 2,000, 3,000, 5,000, 7,000, and 10,000 
cLHS selected samples on 1999/10/31 for the area. It is proved that simulation based on more samples 
can increase the accuracy of SGS simulation, especially under the situation of more than 3,000 samples. 
The simulation with the efficient samples can map the spatial patterns of landscape after Chi-Chi 
earthquake. Figure 8(a) shows the correlation coefficients between the original data and simulation 
data of different sample numbers in seven events. The coefficients of sampling simulation maps are 
close to these of the original data when samples increase. After 5,000 samples, the correlation 
increases slowly. The results indicate that samples number increases a lot in applications, but the 
correlation increasing rate inclines to steady. Figure 8(b) shows the mean absolute error (MAE) 
between the referenced and simulated data. Results are consistent with that indicates the more the 
samples, the higher the match accuracy. SGS simulation based on sufficient samples can capture the 
spatial characteristics of landscape changes including spatial heterogeneity and variability. 

This study also presents the spatial analysis such as the variograms and landscape metrics to explore 
the relationship in different samplings. Figure 9 shows experimental variograms for 100, 500, 1,000, 
2,000, 3,000, 5,000, 7,000 and 10,000 cLHS selected samples and original data on 1996/11/08, 
1999/03/06, 1999/10/31, 2000/11/27, 2001/11/20, 2003/12/17 and 2004/11/19, respectively. These 
experimental variograms show that as the number of samples increased from 100 to 10,000, the ability 
of experimental variograms to capture the spatial structure of actual NDVI images increased. These 
variography results show that the cLHS and SGS approaches can simultaneously select samples from 
multiple NDVI images and simulate all NDVI spatial structures. Results show the variograms are good 
indicators of pattern identification of the NDVI spatial structures in this study. The variation in 
variogram patterns observed among watersheds show that the underlying spatial structure and 
landscape change induced by the disturbances. The variogram using less than 1,000 samples are poor 
spatial representation. From the statistical and spatial analysis, effective samples (5000 samples) can 
be selected by cLHS (only 8% of total) and replicate the statistical distribution and spatial structures of 
the NDVI from the original data. In final, Figure 10 shows NDVI maps produced by SGS simulations 
with 5,000 cLHS samples in the seven events of the area. The SGS results verify that the limits of 
spatial analysis and interpolations of landscape variables are based on semivariograms (or 
autocorrelation functions) solely, stressing the need to account for spatial patterns in highly 
heterogeneous landscapes after large physical disturbances [74]. Therefore, procedures for 
interpolation of NDVI must include information on spatial patterns, either directly from remotely 
sensed images or indirectly by sampling with sufficient spatial variables in the field that depend on the 
interest variable (assuming the sampling is too expensive in the field) [20,74]. The simulated NDVI 
images show that SGS and the cLHS approaches provide effective tools for monitoring, sampling and 
mapping landscape changes. 

To illustrate how to deliver to field surveyors to website, this study set up a prototype based on 
OGC (Open Geospatial Consortium) to display the sampling data distribution. The OGC developed a 
series of standards for geospatial and location based services such as GML (Geography Markup 
Language), WMS (Web Map Service) and WFS (Web Feature Service). 
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Table 4. Statistics of SGS simulations with different samples from NDVI images. 



Samples 


Date 


Mean 


Ma 


Min 


Max 


Skewness 


Kurtosis 


100 


1996/11/08 


0.25 


0.08 


-0.03 


0.44 


-0.40 


-0.42 




1 C\C\C\ /AT lf\/~ 

1999/03/06 


0.24 


0.07 


0.00 


0.41 


-0.42 


-0.27 




1999/10/31 


0.02 


0.04 


-0.12 


0.28 


0.45 


0.52 




r> r\r\r\ / 1 i ir\i-i 

2000/1 1/27 


0.10 


0.04 


-0.06 


0.29 


-0.13 


-0.23 




2001/1 1/20 


0.02 


0.05 


-0.13 


0.46 


1.64 


5.51 




oaat /10/n 

2003/12/17 


A AO 

0.08 


a c\n 

0.07 


A 1 O 

—0.12 


A OA 

0.29 


0.17 


ATT 

-0.33 




OAA/I /I 1 /I Q 

ZUU4/1 


o 7/^ 
U.ZO 


o 

U.UO 


o cxi 
U.U / 


o aq 
U.4y 


n no 
—U.Uo 


n co 

U.J 27 


500 


1 r\r\ s~ it 1 if\c> 

1996/1 1/08 


0.36 


0.05 


0.11 


0.52 


-0.23 


0.07 




1999/03/06 


0.30 


0.04 


0.15 


0.42 


-0.22 


-0.09 




1999/10/31 


0.12 


0.05 


-0.15 


0.26 


-0.42 


-0.14 




2000/11/27 


0.14 


0.05 


-0.05 


0.30 


-0.22 


-0.41 




oaai /1 1 /oa 

2001/1 1/20 


0.29 


0.04 


0.09 


0.46 


-0.15 


0.30 




o a a i / n/ n 

2003/12/17 


A 1 1 

0.13 


A A/I 

0.04 


A AC 

-0.05 


A O 1 

0.31 


A OA 

-0.29 


A 1 1 

-0.1 1 




OAA/I /I 1 /I Q 

ZUU4/1 


o i/^ 
U.Jo 


o 

U.Uj 


A 1 1 

U. 1 1 


o ^7 
U.jZ 


n oi 
— U.ZJ 


n n.7 
—u.u / 


1000 


1996/11/08 


0.35 


0.05 


0.10 


0.49 


-0.83 


1.55 




1999/03/06 


0.26 


0.04 


0.12 


0.41 


-0.04 


0.09 




1 /~\/~\/~\ it /\ /i t 

1999/10/31 


0.13 


0.05 


-0.12 


0.29 


-0.49 


0.09 




2000/1 1/27 


0.13 


0.05 


-0.05 


0.30 


-0.13 


-0.13 




r\ f\f\-\ 1 1 1 /^\ /-\ 

2001/1 1/20 


0.31 


0.06 


0.02 


0.50 


-0.64 


0.61 




2003/12/17 


A 1 O 

0.12 


A A/I 

0.04 


A A/I 

— U.U4 


A OA 

U.29 


A 1 A 

— U.14 


A A/I 

— U.U4 




OAA/I /I 1/1Q 

ZUU4/1 


o ii 
U.J J 


o fK 
U.Uj 


o 1 7 
U. 1Z 


o ^ 1 
U.j 1 


U.Z1 


U.Uj 


2000 


1 r\r\ s~ it 1 if\c> 

1996/1 1/08 


0.37 


0.04 


0.14 


0.52 


-0.56 


0.77 




1999/03/06 


0.29 


0.04 


0.14 


0.42 


-0.1 1 


-0.1 1 




1999/10/31 


0.14 


0.05 


-0.20 


0.31 


-0.70 


0.43 




2000/1 1/27 


0.14 


0.05 


-0.1 


0.35 


-0.10 


-0.36 




r\ f\f\-\ It 1 lr\c\ 

2001/H/20 


0.36 


0.05 


0.07 


0.52 


-0.64 


1.07 




oaat / n/ n 

2003/ 12/ 17 


A 1 1 

0.13 


A A/I 

0.04 


A AC 

-0.05 


ATT 

0.33 


A 1 A 

-0.10 


A O 1 

—0.21 




OAA/I /I 1 /I Q 

ZUU4/1 


U.J / 


o o/^ 
U.Uo 


U. 1U 


U.j4 


— U.ZZ 


— U.ly 


3000 


1 r\r\ s~ it 1 if\c> 

1996/11/08 


0.37 


0.04 


0.18 


0.49 


-0.57 


0.65 




1999/03/06 


0.29 


0.04 


0.15 


0.42 


-0.17 


-0.08 




1 r~v/^/^ it r\ l^> 1 

1999/10/31 


0.14 


0.05 


-0.19 


0.33 


-0.76 


0.58 




^\ /-\ /-\ /-\ / l 1 lr\i—i 

2000/H/27 


0.14 


0.05 


-0.04 


0.32 


-0.18 


-0.29 




2001/1 1/20 


0.37 


0.04 


0.07 


0.51 


-0.59 


0.64 




oaat / n/ n 

2003/12/17 


A 1 /I 

0.14 


A AC 

0.05 


A 1 1 

-0.1 1 


A O 1 

0.31 


A 1 O 


A AA 

-0.09 




O AA/I /I 1 /I A 

ZUU4/1 


o i c 
U.J j 


o oc 
U.Uj 


U. 1Z 


U.jZ 


A 1C 

U. 1 J 


A 1Q 

— U.ly 


5000 


1 r\r\ z' it 1 //ao 

1996/1 1/08 


0.37 


0.04 


0.15 


0.48 


0.58 


0.60 




1999/03/06 


0.30 


0.04 


0.14 


0.43 


0. 1 9 


-0.09 




1999/10/31 


0.14 


0.06 


-0.20 


0.29 


-0.81 


0.61 




2000/1 1/27 


0.15 


0.06 


-0.10 


0.31 


-0.29 


-0.29 




2001/1 1/20 


0.33 


0.04 


0.03 


0.48 


-0.20 


0.16 




TAAl /1 O /1 T 

2003/12/17 


A 1 /I 

0.14 


A AC 

U.Uj 


A 1 A 

— U.1U 


A T 1 

U.31 


a n 
— U.22 


A 1 H 

— U.lo 




70.0-/1 /1 1 /1 Q 

ZUU4/1 


o i c 
U.J j 


o oc 
U.Uj 
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U.1U 


U.jZ 


—U.Zo 


A 1A 

— U.1U 


7000 


1 r\r\ s~ it 1 if\c> 

1996/11/08 


0.37 


0.04 


0.12 


0.49 
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0.81 




1 AAA /AO lf\ S 

1999/03/06 


0.30 


0.04 


0.13 


0.42 


-0.25 


-0.05 




1999/10/31 


0.14 


0.06 


-0.18 


0.32 


-0.84 


0.55 




^\ /-\ /-\ /-\ / l 1 lr\i—i 

2000/1 1/27 


0.14 


0.06 


-0.10 


0.33 
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-0.30 




2001/1 1/20 


0.37 


0.05 


0.04 


0.51 


-0.70 


0.75 




7001/1 9/17 
ZUUJ/ 1Z/ 1 / 


n 1 a 

U. 14 


U.Uj 


— n 1 n 

U. 1U 


n 1 1 

U.J 1 


— n 7A 

U.Z4 


— n H7 
u.u / 




2004/11/19 


0.35 


0.05 
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-0.05 


10000 


1996/11/08 
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0.04 


0.16 
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1999/03/06 
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0.43 
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1999/10/31 


0.14 


0.06 


-0.17 


0.33 


-0.91 


0.77 




2000/11/27 


0.15 


0.06 


-0.10 


0.32 


-0.27 


-0.36 




2001/11/20 


0.37 


0.05 


0.09 


0.50 


-0.75 


0.82 




2003/12/17 


0.14 


0.05 


-0.10 


0.31 


-0.22 


-0.20 




2004/11/19 


0.35 


0.07 


0.06 


0.54 


-0.27 


-0.11 
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Figure 8. (a) Correlation coefficients and (b) MAE between the original data and SGS 
simulation data with different sampling data in seven events. 
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Figure 9. Experimental variograms of NDVI simulations using different samples on (a) 
1996/11/08, (b) 1999/03/06, (c) 1999/10/31, (d) 2000/11/27, (e) 2001/11/20, (f) 2003/12/17, 
and(g) 2004/11/19. 
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Figure 10. SGS simulated NDVI images based on 5,000 samples for area on (a) 
1996/11/08, (b) 1999/03/06, (c) 1999/10/31, (d) 2000/11/27, (e) 2001/11/20, (f) 2003/12/17, 
and(g) 2004/11/19. 
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Figure 11. (a) Sampling data is overlaid on the latest image of Google Earth; (b) The 
server delivered a KML according to end-user's request on Google Earth. 




(b) 

The GML is a geospatial data standard that is neutral to commercial GIS software data format, as 
well as there is more and more GIS software complaining with GML. The WMS and WFS enable data 
providers to publish their by the approach that is based on IT standard and easy to implement. 
Figure 1 1(a) showed the sampling data is overlaid on the latest image of Google Earth. When end user 
sent a WMS request to the server, the server will response the image of sampling data. Google Earth is 
able to overlay the image from WMS response. Thus, end-users can browse the sampling data via 
Google Earth. Moreover, the server can response KML data that is a subset of GML. Figure 11(b) 
illustrated the server delivered a KML according to end-user's request and the KML is overlaid on 
Google Earth. 
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4. Conclusions 

In the paper, multiple NDVI images, which can be generated almost immediately after the remotely 
sensed data are acquired, were used as the detection of landscape changes induced by a large 
disturbance. This study presents an effective framework that integrates cLHS, SGS, and spatial 
analysis in remotely sensed images for efficient monitoring, sampling, and mapping of the impacts 
from chronologically large disturbances on spatial characteristics of landscape changes including 
spatial structure, variability and heterogeneity. 

The cLHS approach is an effective sampling approach for multiple NDVI images from the 
multivariate distributions to replicate the statistical distribution and spatial structures of the NDVI 
images. Using the spatial analysis such as MoranT, and variography, SGS with sufficient samples 
generates multiple realizations and a realization average of NDVI, as well as captures the spatial 
variability and heterogeneity of disturbed landscapes. Spatial analysis of pre- and post-NDVI images 
of a large disturbance is essential for characterizing and quantifying the spatial variability, structure, 
and heterogeneity of landscapes induced by a disturbance. Moreover, landscape metrics have been 
proved effective in land-use change detection because they can characterize the differences between 
the events. The results illustrated that the impacts of large-physical disturbances on spatial variability 
existed and depended on disturbance magnitudes and paths, but were not always evident in 
spatiotemporal variability of landscapes in the study area. 

In sum, the sufficient number of NDVI samples using cLHS (only 8% of total) can be applied to 
monitor the land cover change, which was induced by large physical disturbances. Then, the sampling 
data are demonstrated on the latest image of Google Earth by using Open Geospatial techniques. From 
the spatial analysis, the spatial simulations based on the sufficient samples are found to capture the 
spatial patterns of original NDVI distributions. Land cover changes in a watershed can alter 
hydrological processes and cause severe damages. When the typhoons came in the area, they brought 
about landslides and debris flows. Thus, the detection and concern about land cover change have 
benefit to soil and water conservation. This study will further research on the relationship between land 
cover change and hydrological processes. 
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